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A systematic approach for the construction of a density 
functional for van der Waals interactions that also accounts 
for saturation effects is described, i.e. one that is applicable 
at short distances. A very efficient method to calculate the 
resulting expressions in the case of flat surfaces, a method 
leading to an order reduction in computational complexity, 
is presented. Results for the interaction of two parallel jel- 
lium slabs are shown to agree with those of a recent RPA 
calculation (J.F. Dobson and J. Wang, Phys. Rev. Lett. 82, 
2123 1999). The method is easy to use; its input consists of 
the electron density of the system, and we show that it can 
be successfully approximated by the electron densities of the 
interacting fragments. Results for the surface correlation en- 
ergy of jellium compare very well with those of other studies. 
The correlation-interaction energy between two parallel jellia 
is calculated for all separations d, and substantial saturation 
effects are predicted. 

PACS numbers: 31.15.Ew,71.15.Mb,34.50.Dy,02.70.-c 



I. INTRODUCTION 

The densitypfmictional theory (DFT),§ with its local- 
density (LDA)a'0 ancLsgmilocal generalized-gradient ap- 
proximations (GGA) ,Qi3 is not only successful in numer- 
ous applications to individual molecules and dense solids. 
It is also under intense development, for instance, in order 
to includ&|non-local effects, such as van der Waals (vdW) 
forces. Bilj The latter are needed in order to allow DFT 
to describe sparse matter. A unified treatment of vdW 
forces at .large and asymptotic molecular separations is 
availablCjEJ and a description at short distances and at 
overlap is striven for. An accurate calculalipn for the in- 
teraction of two He atoms has been given,E£l and recently 
the first microscopic (RPA) calculation of the vdW inter- 
action between two self-consistent jellium slabs has been, 
reported and given a density functional (DF) account .til 
The ultimate challenge is to construct an approximate 
vdW DF that is generally applicable, efficient, and accu- 
rate. 

We here propose an explicit form for the vdW DF that 
applies to flat surfaces, test it successfully against these 
slab results, and apply it to two parallel flat semi-inflnite 
metal surfaces. This is a case with relevance for many 
physical situations, including wetting and atomic-force 
microscopy (AFM). Compared to the vdWj-DF approx- 
imation proposed by Dobson and WangjlJ the virtues 
of our functional are the computational simplifications 



gained from choosing a particular sub-class of response 
functions, utilizing a differential formulation and sparse 
matrices, and recognizing the insensitivity to the details 
of the density profiles, simplifications which might trans- 
fer even to three dimensions. 

The ubiquitous van der Waals force plays an important 
role for numerous physical, .cbfrpical, and biological sys- 
tems, such as physisorptioUjESEZl vdW complexes,E3 vdW 
bonds in crystals, liquids, adhesion, and soft condensed 
matter (e.g., bio macr omolecules , biosurfaces, polymers, 
and membranes). liSllj 

The DFT expresses the ground-state energy of an in- 
teracting system in an external potential v{r) as a func- 
tional E[n\ of the particle density n(r),|Which has its min- 
imum at the true ground-state density.EI The Kohn-Sham 
form of the functional makes the scheme a tractable one, 
as it leads to equations of one-electron type, and accounts 
for the intricate interactions among the electrans with an 
exchange-correlation (XC) functional i?xcM-0 This XC- 
energy functional can be cxpxe^cd exactly as an integral 
over a coupling constant (A),l3'Q'EJ the so-called adiabatic 
connection formula (ACF), 
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l)V] - E,,u, (1) 



where V{r, r') = 1/ |r — r'| and where the densibtdensity 
correlation function is denoted by x(r, r', iu; A).l2I EgcU is 
the Coulomb self-energy of all electrons, which is exactly 
cancelled by a corresponding term in x(A, iu)V. Equa- 
tion (U) shows a truly non-local XC interaction and is a 
starting point for approximate treatments, local (LDA), 
semilocal (GGA) and non-local ones. 

The LDA and GGA are completely unable to express 
the vdW interaction in a physically sound way. The exact 
XC energy functional, on the other hand, of course en- 
compasses such interactions .a The basic problem of mak- 
ing DFT a working application tool also for sparse mat- 
ter is to express the truly non-local vdW interactions 
between the electrons in the form of a simple, physi- 
cal, and tractable DF. Equation (|l|) is then the starting 
point. Along these lines, we have proposed extcuidous. pf 
the DFT to include van der Waals interactions,Q'[lliiJlij'E£l 
with very promising results for the interaction between 
two atoms or molecules, an atom and a surface, and 
two parallel surfaces, respectively. This now unified 
approachEj applies for separated systems, i.e. when 
the electrons of the interacting fragments have negligi- 
ble overlaps. 
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The corresponding asymptotic expressions have singu- 
lar behaviors at short sepatajtions d. Yet one knows that 
the vdW forces are finite. They should go smoothly 
over to the XC forces that apply in the interior of each 
electron system. Ihis-phenomenon is often called damp- 
ing, or saturationB^c^ Approximate saturation functions 
have beeiijjroposed, in particular for Jjjia. cases of vdW 
moleculesEil and physisorbed particles 

The key difficulty in extracting the vdW DF from 
Eq. (|^) is the computational complexity. A direct solu- 
tion gives simply too many operations on the computer. 
The guideline for our reduction of the number of such op- 
erations is to exploit analytical advantages of RPA-like 
approximations, to focus on the key quantity, to recast 
the integral formulation into a differential one, leading to 
a sparse-matrix computation, and to make maximal use 
of symmetry. 

Our exploratory study here concerns cases with vdW 
forces between two flat parallel model systems. We first 
test our approximate functional on the model system 
of two self-coasistent jellium slabs, utilizing the recent 
RPA results,L3 which gives the size of the correlation- 
interaction energy per unit area, showing saturation. We 
then test our DF against-aecurate calculations of the sur- 
face correlation energy,r Ir^ showing an excellent agree- 
ment. After these successful tests we make predictions 
on two parallel semi-infinite jellia. 

II. GENERAL THEORY 

The A integration of Eq. ([T]) can be performed an- 
alytically in some cases, such as in the random-phase 
approximation (RPA). In 1957 Gell-Mann and Bruck- 
ner (GMB)c2l presented the RPA correlation energy as 
a selected summation of ring diagrams, which gives a 
logarithmic form.E3 Their study concerns the homoge- 
neous electron gas, where equations simplify thanks to 
the three-dimensional translational invariance and plane 
waves. Here we treat systems with less symmetry. 

By virtue of the fluctuation-dissipation theorem, the 
density-density correlation function x is equal to the den- 
sity change dn induced by an external potential <i>cxt, i-e. 
Sn = x^cxt- It satisfies 

X(A, iu) = xi-iu) + Xx{iu)Vx{\ iu), (2) 

where x is the density response to a fully screened poten- 
tial $, i.e. 6n — x^. We assume here that the coupfing 
dependence of x can be neglected, when performing the 
A integral in Eq. (Q). This is true in the random- phase 
approximation, where x is the density response function 
for A = 0, and is also true for the approximate dielectric 
functions, which we use here. Equation (|l|) then becomes 

E.c = 2^ Re [Tr[log(l-xMV^)]]- itself, (3) 
where the real part means the principal branch. 



To simplify this expression and to get a functional in 
terms of the electron density n(r), we have to focus on the 
key target, the non-local part, introduce key quantities, 
and rewrite the expressions, in order to make physically 
sound and computationally efficient approximations. It 
is more convenient to introduce the polarizability or di- 
electric function instead of x. The polarizability a (a 
matrix in the spatial positions) is defined by the relation 
P = a • E, where P is the polarization. We have 

= -V ■ P = -V ■ a • E = V • a • V$, (4) 

so that from the definition of x, one has x = V • a • V. 
In turn the dielectric function is given by e = 1 + Ana. 
In terms of e Eq. (|) then transforms to 

poo 1 

E.c^ -^Re[Tr[log(V-e-VG)]]-i?,cif, (5) 
Jo ^'^ 

where we have introduced the Coulomb Green's function 
G = -V/A-K and used V^G = 1. The Tr [log] expres- 
sion gives great advantage for the further analytical and 
numerical treatment. The only approximation made so 
far is the neglect of the coupling constant dependence of 
X when doing the coupling constant integration. This is 
not an additional approximation either in the RPA or for 
the approximate e's we use here. 

In order to develop long-range functionals, one may 
substitute approximations for the dielectric function 
based on the free electron gas into Eq. (||). To obtain 
tractable expressions it will normally be necessary to 
make still further approximations. In this case it is de- 
sirable to use the additional approximations only for the 
non-local part of E^^ so as to avoid destroying the ac- 
curacy of the LDA in the high-density regions. Ideally 
one would subtract from Eq. ^ the LDA version of the 
same approximation, and would envisage adding back a 
better version of the LDA. Here we make a similar, but 
more tractable subtraction, which takes the form 

/.OO 7 

El= / r^Re[Tr[log(6)]]-£;,eif. (6) 

E'^^ has the property that it is a good approximation for 
a slowly varying system, becoming exact for a uniform 
system. For density variations slow on the scale of the 
range or width of e(r, r'), it agrees with the LDA, the 
trace in Eq. ^ replacing the integral over density. 
Subtracting Eq. (||) from Eq. (||), one obtains 

/.OO 7 

^"c - i^^^ [Tr [log(e-iV • e ■ VG)]] . (7) 
Jo 

We will call this the non-local exchange-correlation en- 
ergy, although for models more general than those used 
in this paper, an additional short-range correction must 
be applied to make E^^ correspond precisely to the LDA, 
and hence to make E"^' the deviation from the LDA. 

The approximations considered in this paper con- 
tain no non-local exchange component, in effect making 
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Eq. ([7|) our approximation for the non-local correlation 
energy E'^K Using this fact, together with Tr[loga;] = 
log(cleta;) and V^G = 1, we obtain 

/•OC 7 

= X ^log|det(l + 6-i[V,e].VG)| (8) 

where the notation [A, B] means the commutator. We 
will later use the fact that Eq. involves only the de- 
terminant to good advantage. 

III. METHOD FOR PLANAR GEOMETRIES 

Now we are in a situation to discuss what e to use. 
We will in this paper concentrate on the simple case of 
jellium systems. Our aim is to find an ejficient way of 
exploring the planar translational invariance of the jel- 
lium system - not only to decrease the number of spatial 
integrations. The major difficulty in evaluating Eq. (||) 
is the determinant, which is 0{N'^) in the general case, 
N being the number of grid points in a discrete repre- 
sentation. This holds true even in one dimension. So, 
instead of allowing a completely general e, we aim at ap- 
proximations resulting in differential operators only, for 
which the determinant is known to be 0{N), hence sig- 
nificantly simpler to calculate. 

In the particular case of planar translational invari- 
ance, i.e. for planar surfaces or slabs, we use an approx- 
imate form that is made local in the coordinate perpen- 
dicular to the surface, 

e(r,r') = S{z - z') J -0^e^{zy^<^-^'\ (9) 

where k is a wave vector parallel to the surface. Keeping 
the fully non-local form along the symmetry plane al- 
lows, e.g., the effect of the cutoff, which was introdiiecd| 
artificially in previous approximations of this typej^'t^ 
to occur in a natural way laterally. 

The first thing we note about Eq. (^) is that we easily 
form the inverse 

e-\v,v') = 5{z-z') J ^e..(z)-ie^'^-(--'). (10) 

Evaluating the commutator then yields 

.-[V,l = M(.-/i/|5,|M.«'-0, („, 

where the prime indicates differentiation with respect to 
z. In what follows we shall substitute l{z) = log(e(z)), 
yielding l'{z) = e'{z)/e{z). 

In the same basis we express the Green's function, 

" "'^ = / W?^'^' ~ ^')e^'-(--"-'\ (12) 

where 



G,{z-z')^~-e~^V-^-'\. (13) 

Since the logarithm of Eq. (Q) can be expanded in pow- 
ers, the integration over k may be singled out, and we 
may express the non-local correlation energy per surface 
area {A) as 

1"°° rill r rl'^h 

EflA = ^ y _|^log|det(l + /^a,Gfe)|. (14) 

In Eq. (p^, the determinant is given in terms of inte- 
gral operators. To take advantage of the locality of the 
Laplacian, we use {d^ — k^)Gk — 1 to express it in terms 
of differential operators 

det{l + l[d,Gk)^^, (15) 

00 

where 

(I) = det{dl - + l[d^) (16) 

and where 0o is the empty space (e = 1) value of Eq. ([l^). 
The step from Eq. (^J) to Eq. ( p^ ) requires that the 
differential operators are defined throughout the whole 
space. 

Our ultra-fast method is made possible by the obser- 
vation that the determinants in Eq. (|l5|) can be written 
down, not only for the full system, but also for a sub- 
division of it. Related determinants for the subsystem 
satisfy a simple second-order differential equation as a 
function of subsystem size. Thus by a simple renormal- 
ization, one may evaluate i?"' with the same effort as 
finding the charge induced by an applied electric field. A 
similar relation holds also in several dimensions, which 
will be explored in another paper. 

To make this more concrete, let us suppose that efe(z) 
varies only in the interval < z < L (which will eventu- 
ally be extended to infinity) and takes the same value at 
either end point. This is the case for parallel surfaces or 
slabs of identical materials. Then for each value of z we 
can define a determinant 4'{z) for the subsystem extend- 
ing from to z. It is clear then from Eqs. (|l^, ([T^), and 
(O) that is given by 

ni / . , f°°du f d^k , <p(L) , , 

£;,"'M= hm / — / log-p-i. 17 

l^coJq 27r J {2ny 4>o(L) 

As discussed in Appendix A, the determinants (j){z) and 
00 (^) individually have oscillating signs that do not occur 
in their quotient. However, the envelope determinants 
0(z) and (j>o{z) can be scaled so that they satisfy the 
simple differential equation 

{ek4>'y = khk4>, (18) 

together with the boundary conditions that 0(0) = and 
0(L) = 1. In terms of 0, we obtain 
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Ef/A=-iu^ 



d^k (j)'{0) 



(19) 



where the prime indicates differentiation with respect to 
z, which in this case is the subsystem size. However, 
note that </) is also just the electrostatic potential as a 
function of distance z, within a system having a potential 
difference across it along with a prescribed variation in 
e. Thus the calculation of the determinant becomes a 
simple electrostatic problem which is easily solved. 

To illustrate this, consider the case of two identical par- 
allel surfaces a distance d apart, when d is much larger 
than the thicknesses of the surface-healing layers. Solv- 
ing Eq. (^8|) for the described boundary conditions then 
becomes a trivial matching problem (see Appendix B), 
which after insertiopj|Bto Eq. (^|) immediately leads to 
the Lifshitz formulMEl 



^ni,, rdu f d^k ^ . 



27„i. (20) 



Here p — {cb — l)/{€b + I), eb being the bulk dielectric 
function, and 7ni is defined by 



7„i = (E^^'id ^ ^) - E^\0))/2A 



(21) 



Since, by construction, i?"' = for a uniform {d — 0) 
system, 7^1 may equivalently be defined as the non-local 
correlation contribution to the surface tension of a single 
surface. 

The original ACF Eq. (|l]) is now reduced to a set of 
simple electrostatic calculations, each one being an 0{N) 
operation instead of 0{N'^), a major simplification. Of 
course the success of the method depends on how well we 
can reproduce the true dielectric function using our ap- 
proximate form Eq. (^. The only approximation made 
so far is the assumption of a local dielectric function per- 
pendicular to the surface. 



IV. APPROXIMATE DIELECTRIC FUNCTION 

Equation (|l^) or ( pj| ) provides the basis for a func- 
tional that describes the van der Waals interaction be- 
tween planar objects. To turn these equations into den- 
sity functionals, we have to introduce quantities that de- 
pend on the density n(r). Our suggestion is based on an 
approximate dielectric function ek that depends on the 
local density n{r). It utilizes experiences from the ho- 
mogeneous electron gas and from experimental studies of 
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e(q,uj) 



the dynamical structure factor S'(q, uj) oc Im 

where ?iq and huj are the momentum and energy losses, 
respectively, of a photon or a charged particle being scat- 
tered while passing a bulk sample. There is a peak in 
^(q, w), the plasmon peak, sharp in the ideal electron 
gas and of varying width in real materials. This peak 
carries most of the spectral strength and has lo equal to 



the plasma frequency as g — > 0, and then a disper- 
sion with a limiting behavior ujq — > /2m, the kinetic 
energy of one electron, in the impulse-approximation, 
valid in the Compton-scattering limit, q —f 00. In the 
electron-plasmon coupling one focuses on the magnitude 
and position of the sharp plasmon peak, and neglects 
the broadeningr-ii-e. e is described in a plasmon-pole 
approximation.L3 A dispersion law like 



(22) 



has been shown to efficiently account for the average 
behavior of plasmon-like excitations and foii-correlation 
properties of the homogeneous electron gaso Introduc- 
ing the electron density via u>p{z) = Airn^z) and vp(z) = 
(377^71(2))^/"^ in Hartrec units, n{z) being the electron 
density profile in this planar case, the dielectric function 
can be written as 



e(z, fc, iu) = 1-1- 



uj^iz) 



u2 + {vF{z)q{k))y3 + q\k)/4' 



(23) 



where the imaginary frequency w = iu is introduced. 

Alternatively, Eq. ( p3| ) may be viewed as an interpola- 
tion between the exact small- and large-q behavior of the 
Lindhard expression for the frequency-dependent dielec- 
tric function. However, we see little point in using such 
an elaborate expression, since our concern here is to in- 
vestigate how well local approximations to the dielectric 
function work in highly non-uniform systems. 

In directions parallel to the surface, our approxima- 
tion Eq. ( p3| ) allows fully for the non-diagonality of e 
with respect to the corresponding spatial coordinates, as 
implied by a Fourier transform with respect to the par- 
allel wave vector k. However, in directions perpendicular 
to the surface, our approximation takes e to be diagonal 
in the coordinates z and z' . It is thus taken to be local 
not only in this sense, but in the additional sense that it 
is a function only of the local density. To compensate, 
we retain the corresponding component q± of the wave 
vector q in the right side of Eq. (|2^) as a parameter, so 
that everywhere = fc^ + . We thus take l/q± to be 
a constant measure of length over which e is effectively 
nonlocal. The dispersion perpendicular to the surface in 
Eq. (|2^) is in this way replaced by a parameter that we 
will fix to some length scale appropriate for the surface. 

For physical reasons such a length scale should be as- 
sociated with intrinsic electron-gas parameters like the 
screening length or the extent of the correlation hole. 
There are of course several choices available. It must be 
kept in mind that we are after long range surface prop- 
erties in a variety of environments. These properties are 
determinedJpy various response functions introduced by 
FeibelmanJ^ of which the simplest. 



d(iu) = 



Jdz znind(iu, z) 
Jdz Ui^diiu, z) ' 



(24) 



is the centroid of induced charge when a uniform electric 
field is applied perpendicularly to the surface. However, 
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for the van der Waals properties of a planar s tir race, a 
related function D{iu) as defined by Hult et aZ.,Eil 



D{iu) 



eb{iu) + 1 efc(m) + 1 



d{iu), 



(25) 



is more important .E3E3t£3'liJ In particular the D-function 
arises in connection with the calculation of van der Waals 
planes, which are determined not only by the D of the 
surface in question, but also by a response function of 
the other body. For example, for a surface in the vicinity 
of i j^ oj ^pic atom, the van der Waals plane Z is given 



1 



^vdW 



dua{iu)D{iu), 



(26) 



where a{iu) is the polarizability of the atom and C3 is 
the coefficient of the leading term in the asymptotic 
form of the interaction energy. In order for the surface 
calculation to scale correctly for a wide variety of atoms 
with different Q;(iu)'s, it is obviously important for our 
approximation to well reproduce D{iu). Similarly, for 
two parallel surfaces labeled A and B, the van der Waals 
plane for surface A is given by 



vdW 



(47r)2C2 



dup^{iu)D^(iu) + AZ^, (27) 



where p^{iu) is the long wavelength surface response 
function of surface B, p-^(m) = (ef (m) — l)/(ef (iu) + 1), 
and C2 is the coefficient of the leading 1/ term in the 
asymptotic form of the interaction energy. The main dif- 
ference in cases such as this one, where two infinite bod- 
ies are involved, is that terms involving multiple reflec- 
tions no longer vanish in the asymptotic limit, and should 
be explicitly included as indicated by the AZ"^ term in 
Eq. (|27|). However, it is common experience that in the 
asymptotic limit multiple reflection terms arp. .typically 
small enough to treat in perturbation theory,E3Ej a con- 
clusion that we agree with. Therefore, the first term in 
Eq. (p7| ) is dominant, and the conclusion that D{iu) is the 
key quantity to obtain accurately remains. Note, how- 
ever, that in the numerical calculations presented later, 
the full form of Eq. ( p7| ) , including all multiple refiections 
(see Eqs. (8) and (9) of Ref. is used. 

Thus we opt to choose q± based on the premise that 
the D-function should be reproduced accurately. We im- 
plement this by choosing q± so that D{0) agrees exactly 
with full LDA calculations of this quantity, a procedure 
precisely analogous to that of Ref. |l^. In this way, the pa- 
rameter q± may be indirectly determined as a function of 
the electron gas parameter r^. This determination gives 
q± ^ 1/rs ^ kp as expected from previous arguments. 

The constant q± smoothly limits the response at small 
k and small u values. It thus-xeplaces the sharp cut- 
off used in the earlier schcmcEj In doing this, we con- 
sciously violate the Lifshitz limit, obtaining somewhat 
smaller vdW coefficients than Ref. O, since the choice 



q± > affects the fc ^ limit of Eq. (g3|). This also 
means according to Eqs. (26) and (|27| ) that the van der 
Waals planes will be predicted to be somewhat too far 
from the surfaces in this approximation. The d-function, 
i.e. Eq. ( p^ will also be too large, significantly so at 
small u. 

Far more important, though, is the fact that our sim- 
ple dielectric function reproduces the dynamic properties 
of the £)-function very well, as shown in Figures 0and|. 
Thus the overall scaling properties of our theory in a va- 
riety of non-uniform van der Waals environments should 
continue to be correct. 
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FIG. 1. The D- function D(iu) for a jellium profile of 
Ts = 2. Solid line: Present calculation. Dotted line: 
Reference W Circles: D-function based on a TDLDA 
calculationita 
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FIG. 2. The _D-function D(iu) for a jellium profile of 
Ts = 4. Solid line: Present calculation. Dotted line: 
Reference W Circles: D-function based on a TDLDA 
calculation.tJ 
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FIG. 3. Optimal values of q±, calculated for jellium sur- 
faces of different densities, described by Vg. Circies: Optimal 
q± as a function of r^, calculated from LDA dataCj Solid line: 
Equation 

In order to investigate the approximation with regard 
to key properties, we use as input to our functional a 
set of self-consistent single-surface jellium densities in the 
metallic range (r^ — 2 — 5). The exact form of the density 
profiles is however not that important; by construction, 
a non-local energy is expected to depend less sensitively 
on the exact form of the density, and results in this paper 
show that this is indeed the case. 

The optimal qj_ values found by fitting to the values 
of D{0) of Ref. II are shown in Fig. |^ as a function of 
the electron-gas parameter r^, well accounted for by the 
simple interpolation formula 

qi_ = 0.416e-"-2i7'-= + 0.168. (28) 

The variation is roughly 50 percent over the whole metal- 
lic range, indicating a rather small overall sensitivity. 

Table | shows how the quantities D{0) and Z^dw 
varies with r^, using the parameterization Eq. (p8|). In 
the fourth and fifth columns, the values of the present 
method for the van der Waals plane Z^^w are compared 
to those of Ref. |l^. Although the slightly larger coef- 
ficients are expected from considerations earlier in this 
section, the fact that the numbers come out almost in- 
dependent of Vs is somewhat of a surprise. 



V. RESULTS 

To test our approximate DF, we solve Eq. ([l8|) , us- 
ing Eq. (|2^), and insert the result into Eq. (|19|) for a 
known system, consisting of two parallel jellium slabs, 
separated by a distance d. Figure shows results in 
ergs/cm^ (ergs/cm^ ~ 0.6423 /xHa/a^) for the non-local 
correlation-interaction energy per surface area, {E^^{d) — 
E:^^{oo))/A. The results using the — 2.07 value of 



TABLE I. Static and dynamic image-plane positions (a.u.) 
of the jellium surface, calculated as a function of r^. 





D(0)" 


D(O)'^ 




2.00 


1.57 


1.57 


1.15 


2.07 


1.55 




1.15 


2.30 


1.48 




1.15 


2.66 


1.41 




1.14 


3.00 


1.35 


1.35 


1.13 


3.28 


1.32 




1.14 


4.00 


1.24 


1.25 


1.12 


5.00 


1.17 


1.17 


1.09 



0.96 

0.87 
0.79 



^This work. ''Reference |o[ '^This work. '^Reference |l|. 

Eq. (^8|) are compared to those of a recent RPA cal- 
culatioiii of a slab system of Ts = 2.07 by Dobson and 
Wang,til as well as to those of their approximate DF 
(IGADEL). The saturation effects are found to be sub- 
stantial in this small-separation region, and we judge all 
the proposed density functionals in the figure to give good 
accounts of the non-local correlation energy. The agree- 
ment with IGADEL reflects the inherent similarities be- 
tween IGADEL and our approximation. The tractability 
of the latter is however not reflected in the table, but 
has to be stated here (about a thousand times faster 
than IGADEL for this particular system, due to the 
overall lower computational complexity of our DF), to- 
gether with the claim that this gives great prospects for 
a tractable future general DF. 

The calculation presented in Fig. ^ is performed using 
a simple superposition of the densities of two separate 
slabs (obtained from d = 12 data) as input. The dif- 
ference between the non-local correlation energy for the 
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FIG. 4. Small-separation (d) variation of the non-local 
correlation- interaction energy (ergs/cm'^) betwepi two paral- 
lel jellium slabs of Ts — 2.07 and width 5 a.uO Solid line: 
Present calculation. Stars: Same calculation but using the 
self-consistent densities of Ref. ^[ Circles: RPA— LDA of 
Ref. |l]. Diamonds: IGADEL-LDA of Ref. |l|. 
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TABLE II. Non-local correlation contribution to the sur- 
face energy of the jellium surface, 7ni (ergs/crn^, as func- 
tion of Ts ^calculated from Eqs. ( |l9| ) and (^3|) together 
with Eq. (ESI) and self-consistent single-surface LDA densi- 
ties. Comparison is made with results for 7ni obtained by 
subtra|e*jng the LDA contribution to the surface correlation 
energyEII from reftwlts of self-consistent RPA calculations for 
a jellium surface.Ea 















2.00 


451 


476 


2.07 


408 


435 


2.30 


300 


332 


2.66 


196 


226 


3.00 


138 


163 


3.28 


106 


129 


4.00 


60 


74 


5.00 


31 


39 



^Non-local correlation energy for single-surface jellium. 
'^RPA-LDA (in the RPA) of Refs. || and 



self-consistent density of the slab system and that ob- 
tained by superposition turns out to be very small, as 
indicated by the stars in Fig. ^. 

The minor deviations at d close to zero between our 
results (solid curve) and those of the full RPA calculation 
(circles) are due to our somewhat crude treatment of the 
width of the exchange-correlation hole perpendicular to 
the surface, a small price one has to pay for a tractable 
DF. 

Another important property is the surface correlation 
energy. This quantity has recently-been calculated for 
a jellium surface within the RPAjHl an approximation 
thought to give the long range correlation effects accu- 
rately. The long-range part (7»j4 may be extracted using 
the data of Kurth and Perdew,EZl by subtracting the LDA 
contribution to the same quantity. In Table ||, the result 
(column 3) is compared with our approximation (column 
2), calculated as the non-local correlation energy given by 
Eqs. (^ and (p3|), for single-surface jellium at various 
values of r^, using Eq. (p8|). We note that the two ap- 
proximations differ on average by only 13 percent. This is 
perhaps the strongest indication that indeed the physical 
approximations made in this paper are robust. 

A remarkable fact already indicated in Fig. ^ is that 
the non-local correlation energy is quite insensitive to 
the exact form of the density profile. To further test this 
assumption, we use a linear superposition of two iden- 
tical self-consistent LDA single-surface densities (SLDA) 
to calculate the non-local correlation surface energy ac- 
cording to Eq. (|2l|). The result closely follows the column 
2 result of Table ||, with a mean error of only 3 per- 
cent. This observation adds to the accumulated findings 
supporting the use of superpositions of single-fragment 
densities in a future general DF. 

After these successful tests of the predictive power of 
the DF, defined by Eqs. (|l|), (|2|), and (H), appHca- 
tions to other systems where no other results are avail- 
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FIG. 5. Small-separation (d) variation of the non-local 
correlation- interaction energy (ergs/cm^) for two semi-infinite 
jellium surfaces of = 2.07 (SLDA density; solid curve), cal- 
culated from our DF defined by Eqs. (^, (^, and (Q. 
The results are compared to the corresponding interaction 
between two 5 a.u. slabs with edge separation d (same as in 
Fig. ^; dashed curve). The dotted curve shows the asymp- 
totic form of the interaction, E — ~C2/(d — 2Zvdw)^, with 
C2 = 1.34 fiRa. and Zvdw = 1-15 a.u., calculated as in Ref. 
but with the dielectric function Eq. (p^. 

able might be done with confidence. Here we present 
results of an application to two semi-infinite jellia of iden- 
tical Tg, having their parallel surfaces a distance d apart 
(Figs. H and g). The input densities are obtained by 
linear superposition of two self-consistent single-surface 
LDA densities (SLDA). We stress that our DF is very 
fast; obtaining a single value for a given density profile 
and a given separation only takes a few seconds on a 
typical workstation of today. 

Figure ^ shows the variation of the calculated non- 
local correlation-interaction energy for two semi-infinite 
jellia of Ts = 2.07 as a function of the separation d. It 
illustrates two facts: the significant deviation from the 
corresponding quantity for two thin jellium slabs (same 
as in Fig. ^; 5 a.u. wide) and the substantial saturation 
effects, the latter by comparing with the results of the 
asymptotic DF formula,li3 applying the dielectric func- 
tion Eq. (|2|). 

In Fig. o, we present the normalized non-local 
correlation-interaction energy {E{d) — E{oo))/2Ajn\ h^ 
the metallic range, showing how the interaction varies 
with Ts. The curves depend only weakly on when 
scaled in this manner. 



VI. CONCLUDING REMARKS 

In summary, we have studied the basis for a DF ac- 
counting for vdW interactions, by starting in the manner 
of our previous work but with essential generalizations 
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FIG. 6. Small-separation (d) variation of the normalized 
non-local correlation- interaction 

energy, {E^^{d) — i?"'(oo))/2A7ni, for a range of values. 
The 7ni values are presented in Table |l[ 



to small separations between interacting objects. A sys- 
tematic approach for the construction of such a DF is 
described, together with a very efficient method to calcu- 
late the resulting expressions. In the case of flat surfaces, 
results for the interaction of two parallel jellium slabs axe 
shown to agree with those of a recent RPA calculation,Eil 
and we show that input densities can be successfully ap- 
proximated by a superposition of the electron densities of 
the interacting fragments. Results for the surface energy, 
of jellium are compared favorably with other studies.EZI 
As a prediction of the theory, the interaction energy be- 
tween two parallel jellia is calculated for all separations d 
and in the whole metallic range. The well known asymp- 
totic behavior {E oc 1/z^) is obtained for large d, and 
as d becomes smaller, substantial saturation effects are 
predicted. 

The major significance of these results is the demon- 
stration that such numbers can be calculated accurately 
at a reduced computational complexity and hence greatly 
improved speed. We have shown that for a subclass 
of dielectric functions, the resulting expressions for the 
non-local correlation energy may be calculated very ef- 
ficiently, and that even a simple approximation to the 
dielectric function yields valuable insight and reproduces 
several physical properties of flat surface and slab mod- 
els. Furthermore, we have indicated that generalizations 
to three-dimensional systems are possible, and that the 
results here suggest such an attempt to be a fruitful 
one. In this way there should be a basis for applica- 
tions to numerous physical, chemical, and biological sys- 
tems, such as vdW bonds in crystals, liquids, adhesion, 
soft condensed matter (e.g., biomacromolecules, biosur- 
faces, polymers, and membranes), and scanning-force mi- 
croscopy. 



ACKNOWLEDGMENTS 

We thank J. Dobson for providing us with several nu- 
merical density profiles obtained in Ref. ^ for parallel 
slabs, and J. Perdew for providing computer code to gen- 
erate density profiles for isolated jellium surfaces. Work 
at Rutgers supported in part by NSF Grant DMR 97- 
08499. Financial support from the Swedish Natural Sci- 
ence Research Council and the Swedish Foundation for 
Strategical Research through Materials Consortium no. 
9 is also acknowledged. 



APPENDIX A: DETAILS OF THE EVALUATION 
OF THE DETERMINANTS </> AND 0o 

Here we show how ratios of determinants like Eq. (|lj) 
can be efficiently evaluated. To eliminate the oscillating 
sign problem mentioned in the main text, we need to first 
go to a discrete representation for the operator {d^ — 
k"^ + I'l^dz) occurring in Eq. Specifically we take 

N points between and L (for the full determinant), n 
points between and z (for subsystem determinants), 
so that z = {n/N)L, with the spacing between points 
h — L/N. We use a similar representation for the empty 
space operator (9f — fc^). Thus, replacing the subsystem 
determinant ^(2) by the discretisized we have 



det 



/ ai 61 
ci 02 

■■• '■. 



\ 







\0 ... c„„ 



(Al) 



In / 



where 



K = l + ^l'k,n 



2 



(A2) 



along with a similar relation for the empty space deter- 
minant 00,71 ■ Being tridiagonal, the determinant can be 
evaluated in 0{N) time, as contrasted to the 0{N^) time 
applicable for a general determinant. 

The determinant is readily expanded by minors, giving 
the recursion formula 



■1+1 — CLn+l4>n — bnCnCj)n~ 



(A3) 



with <f)Q — 1 and (pi = ai. It is clear from (A3) and 
(|A2| ) that (f>n oscillates in sign from order to order, an 
oscillation that is cancelled in the ratio Eq. (15) by a 
similar oscillation in 4>q. 

However, the envelope determinant (— 1)"0„ satisfies 
a simple differential equation in the continuum limit. To 
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make the exact form of this differential equation identical 
to Eq. (^Sj), we further scale the envelope determinant 0„ 
through 

0„ = (-1)"5„0„, (A4) 
where the scaling function Sn = Hf^ibi satisfies 

Sn = bnSn-1 (A5) 



with ^0 = 1. Using (A2), we see that in the continuum 
/i — » limit, (fA^) becomes 



(A6) 



with the boundary condition ^(O) — 1. Equation ( [A^ ) 
has the solution 



S{z) 



(A7) 



This means that if efc(L) = efc(O), then the scaling due to 
S has no effect on the final result. Otherwise, the argu- 
ment of the logarithm in Eq. (^9|) should be multiplied 
by S{L) obtained from (|A7|). 



Th e di ffere nce equation for is obtained by use of ( A4) 
and (A5) in (A3), yielding 



->n+lVn+l 



+ fln+l^ri + Cn<j)n~l = 0, 



with 



and 



ho = -(l + ai/&i). 



(A8) 
(A9) 
(AlO) 



Substitution of ( |A2D into (|A^) yields 

((^n+l - 20„ + (hn-l) - h'^k'^4>n 
+ ^^fe,n+l('^n+l - <i^n-l) = 0. 

In the contimmm [h —>■ 0) limit, this becomes 



(All) 



(A12) 



This equation is the same as Eq. ([T^) in the main text, 
although there we used the notation for a particular 
solution, while in this appendix it represents the general 
solution. 

This general solution to Eq. (|lj) can be written in the 
form 



(z) = a4>a{z) + 134) piz) 



(A13) 



where (t>a{^) = 1 and <j>a{L) — 0, while <j>p{0) = and 
<j>p{L) = 1. Th e bo undary condition (|A9| ) implies that 
a = 1, while ( A10| ) combined with ( |A2D implies that 



4>'{0) = l/h as h approaches zero. This means that 
can be neglected, since (3 will be of order l/h. Specifically 
we have ^'(0) ^ P(t)'piO) = l/h, so that 



(AM) 



The l arge coefficient (oc l/h) cancels out in the ratio of 
( A14) and the analogous free-space expression, so we may 
write 



ML) ML) 



S{L) = 



(0) 



S{L), 



(A15) 



where the continuum version of (A4) was used to obtain 
the first identity, after noting that the oscillating signs of 
the numerator and denominator cancel in the ratio, and 
that the free space value of S is unity. F ina lly substitut- 
ing this into Eq. ([l7|), we obtain, using (A7), 



E. 



nl 



' du 
2^ 



■logT 



^'(0)7^ 



(2^)2 ^4>'(Q)^/7m' 



(A16) 



This is our most general result, which reduces to Eq. (|l9| ) 
when efc(O) = tk{L). Note that in the main text, we used 
the notation (p and (pQ for the particular solutions that 
are called 4>fj and (pa^fj here. 



APPENDIX B: DETAILS ON THE LIFSHITZ 
LIMIT 

Let efe(z) = 1 in between two surfaces a distance d 
apart, and let ek{z) = tb in the bulk of each surface. 
Furthermore, let d be large enough so we can assume 
sharp boundaries between the three regions. Then, the 
interaction energy (£'"'(d) — i?"'(cx)))/A may be calcu- 
lated exactly. Note that although the interaction energy 
may be calculated this way, the constant component con- 
tributing to the surface energy may not, but needs a much 
more involved calculation. 

The matching problem becomes solving Eq. ( p^ ) for cp 
and (po in the regions left bulk (lb), middle region (m), 
and right bulk (rb). Let cpQ = e*^^, with the origin lying 
on the boundary between the left and middle regions. 
Now we want to find the (p that approaches zero far into 
the left bulk, and e^^ far into the right bulk. In the left 
bulk, we must have 



ae 



kx 



(Bl) 



a being an arbitrary constant. Note that the wanted 
quantity (/)g(O)/0'(O) equals 1/a. On the boundary be- 
tween the left bulk and the middle region, (p must be 
continuous, and (p' must have a discontinuous step with 
the size of reflecting the fact that the displacement 
field eb<p' must also be continuous. The solution in the 
middle region now becomes 
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a(cosh(fcz) + ei, sinh(A:z)). 



(B2) 



The same matching conditions apply between the middle 
region and the right bulk, although here, we only need 
to consider the solution growing to the right, since we 
want to match (p to the value of (po infinitely far into the 
right bulk. Obeying the matching conditions means the 
solution in the right bulk becomes 



,growm{ 
rb 



^ (j)n,{d)k + (j)[^{d)/et 
2k 



(B3) 



Equating 00 with (B3) determines the coefficient a, and 
the solution becomes 



</>UO)/0'(O) ^ 



1 



2 

p e 



-2kd 



2kae 



kd 



1 



(B4) 



with p = {e}y — l)/(eb + 1). It is clear from (B4) that the 
constant contribution to log((/)o(O)/0'(O)) and hence to 
the surface energy is given by — log(l — p^). As discussed 
earlier in this appendix, that is not the correct constant 
and should be excluded, yielding Eq. (pO[). 
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